Synchronization of delayed coupled neurons in presence of inhomogeneity 
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In principle, while coupled limit cycle oscillators can overcome mismatch in intrinsic rates and 
match their frequencies, but zero phase lag synchronization is just achievable in the limit of zero 
mismatch, i.e., with identical oscillators. Delay in communication, on the other hand, can exert 
phase shift in the activity of the coupled oscillators. In this study, we address the question of how 
phase locked, and in particular zero phase lag synchronization, can be achieved for a heterogeneous 
system of two delayed coupled neurons. We have analytically studied the possibility of inphase 
synchronization and near inphase synchronization when the neurons are not identical or the con- 
nections are not exactly symmetric. We have shown that while any single source of inhomogeneity 
can violate isochronous synchrony, multiple sources of inhomogeneity can compensate for each other 
and maintain synchrony. Numeric studies on biologically plausible models also support the analytic 
results. 



I. INTRODUCTION 

Synchronous firing of neurons has received much atten- 
tion in relation to the generation of brain wave rhythms 
and information processing at various aspects in the neu- 
ronal systems, such as selective attention and the bind- 
ing problemii^. Although synchronization with nonzero 
phase lag has been observed between distal areas in hu- 
man brains^, zero phase lag synchronization is of par- 
ticular interest with respect to the binding hypothesis'*. 
The mechanism of these phenomena has been subject 
of controversial debate in a more general context; be- 
yond its functional relevance, the zero time lag synchrony 
among such distant neuronal ensembles must be estab- 
lished by mechanisms that are able to compensate for the 
delays involved in the neuronal communication. Laten- 
cies in conducting nerve impulses down axonal processes 
can amount to delays of several tens of milliseconds be- 
tween the generation of a spike in a presynaptic cell and 
the elicitation of a postsynaptic potential^. The ques- 
tion is how, despite such temporal delays, the reciprocal 
interactions between two brain regions can lead to the 
associated neural populations to fire in unison. 

It is quite probable that a variety of mechanisms are 
responsible for bringing synchrony at different levels (dis- 
tinguishing for example, among local and long-distance 
synchrony) and different cerebral structures. Yet there 
are strong evidences that long distance inter-hemispheric 
phase locking at zero phase lag arises via reciprocal con- 
nectivity rather than via locking to a common inpuli^. 
Connections between different areas and cortical columns 
are mediated by excitatory interactions. It has been ar- 
gued that delayed excitatory connections between neuron 
models do not readily synchronize the neurons and in fact 
they usually lead to antiphase firing->S-. This belongs to 
the mechanism by which model cortical neurons make the 
transition from excitable state to repetitive firing. Phase 
resetting curves (PRCs) keep track of how much an input 
advances or delays the next spike in an oscillatory neuron 
depending upon where in the cycle the input is applied. 
PRCs are formally found by perturbing the oscillation 



with a brief depolarizing stimulus at different times in 
its cycle and measuring the resulting phase shift from 
the unperturbed system2r2i. There was an agreement 
that the type of bifurcation which results in repetitive 
firing of the neuron, determines the response of the neu- 
ron to the external stimulations, so the functional form 
of PRC: In type I neurons with mainly positive PRC, 
transition to repetitive firing occurs through a saddle- 
node bifurcation on invariant circle (SNIC) and for type 
II neurons which have PRCs with a large negative lobe, 
repetitive firing ensues birth of a limit cycle through a 
Hopf bifurcatioo^ dPi^^ . Just recently it is shown that 
type II PRCs can also occur in systems that are arbitrar- 
ily close to a SNIC bifurcatior^i^. 

In general the form of the interaction between oscilla- 
tors together with their intrinsic response (the PRC) pro- 
vide sufficient information about the ability of the system 
to synchronize or desynchronize the oscillationai^ii^. On 
the neurons side, the type II oscillators can more eas- 
ily synchronize: A PRC which contains both negative 
and positive lobes can allow inputs to both slow down 
the oscillator which is ahead and speed up the oscilla- 
tor which is behind. Indeed it turned out inphase firing 
of type I neurons with excitatory synapses is only pos- 
sible with unrealistic instantaneous couplings^ii^"— . In- 
tuitively, phase lags equal to the conduction delay are 
expected in causal limit when the delayed arrival of an 
input triggers a spikei^. On the effect of the type of in- 
teraction, it is believed that inhibition and not excitation 
is responsible for synchronization of the neurons in many 
brain areas^. While instantaneous inhibitory couplings 
lead to antiphase evolution of relaxation oscillators, in 
presence of delay inphase firing of is more likely to occur 
with inhibitory synapse o'^'^'-'^^d^ . 

Synchronization is also affected by the configuration 
through which the neurons interact. A prominent exam- 
ple where a special configuration of the links is exploited 
to synchronize delayed coupled oscillators, is suggested 
by Fischer et al.^°. In this arrangement two distant oscil- 
lators communicate indirectly through a relay oscillator 
and the system exhibits zero-lag synchronization between 
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two outer elements provided the two branches are simi- 
lar. While the symmetry is held (two branches are sim- 
ilar) the model is generic and synchronization is almost 
independent of the intrinsic parameters of the outer and 
the relay neurons and also the type of interaction. Yet, 
as the authors indicate the results critically depend on 
the symmetry of the parameters of lateral neurons and 
connections22i2i. 

Presence of inhomogeneity in the networks of neural 
oscillators can also destabilize synchronyi^j^. In classi- 
cal models for synchronization coupling strength and in- 
homogeneity are competing factors which determine the 
collective dynamics of the systemSlsM^ jjj f^ct coupled 
oscillators can overcome finite mismatch in intrinsic fre- 
quencies and match their frequencies, but the cost is a 
finite phase lag^**. How would the nonidentical neurons 
synchronize when they communicate from distant? We 
follow this question in a more general scheme by study of 
two model neurons coupled through delayed connections, 
when the parameters of either the neurons or the connec- 
tions are slightly different. Using iterative maps we give 
a general analytic framework for the existence and sta- 
bility of phase locked solutions for two phase oscillators 
connected by delayed pulsatile couplings, in presence of 
inhomogeneity. The results coincide with the those of 
other authors in special cases, e.g., for two similar os- 
cillators coupled via symmetric reciprocal connections^^. 
The numerical results given both for the phase model and 
for more biologically plausible conductance based models 
for neurons, support the the analytic findings. 

II. MODELS AND METHODS 

To model the type I and type II neurons, we have 
considered Wang-Buzsaki (WB)^^ and classical Hodgkin- 
Huxley (HH)^^ models, respectively. In type I neurons 
with mainly positive PRC, transition to repetitive firing 
occurs through a infinite period bifurcation and with a 
slowly increasing applied current, the neuronal dynamic 
changes from stationary to oscillatory with arbitrarily 
small frequency. For type II neurons which have PRCs 
with a large negative lobe, repetitive firing occurs via 
birth of a limit cycle through a Hopf bifurcation, and 
the onset of the repetitive firing occurs with non-zero 
frequency^ (see also Ref. [13]). 

Both the WB and HH models when stimulated with 
a supra-threshold constant input current settle into a 
limit cycle, embedded in its 3- and 4-dimensional phase 
spaces, respectively. The current balance equation for 
both model neurons is 

= -gnam^Hvi - Vna) - 9kn^{Vi - Vfe) 

at 

-gi{v,-Vi)-hj+h, (2.1) 

where i and j are chosen from (1,2). C is the membrane 
capacitance in ^F/cm?, Vi is the membrane voltage in 
mV , and li is the (density of) applied current in ^A/ cm? . 



The parameters gnai Qk and gi are the maximum conduc- 
tances per surface unit for the sodium, potassium and 
leak currents and Vna, Vk and Vi are the corresponding 
reversal potentials. All the constant parameters for both 
neuronal models are given in TABLE. H] Details of the 
two models can be found in Appendix. 

TABLE I: The parameters for Wang-Buzsaki and Hudgkin- 
Huxley neurons 
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The neurons are assumed to communicate through 
chemical synapses modeled by 

= 9ijSij{t - nj){v, - Esyn), (2.2) 

where g^j is the synaptic maximum conductivity and Ty 
is the delay in communication of the neurons j (presy- 
naptic) and i (postsynaptic). Synaptic reversal poten- 
tial Esyn determines the excitatory /inhibitory type of the 
synapse. Throughout this manuscript, we take E^yn — 
lOmV for excitatory and Egyn — 70mV for inhibitory 
synapses, respectively. The rate of change of the synap- 
tic variable s(t) is given by the following equation 

= af{vj - vth){l - s^j) - ^Sij, (2.3) 

with a and /3 defining the synaptic activation and deac- 
tivation time constants, respectively, and f{x) = 1/2 [1 -|- 
tanhirjx)] ensures activation of the synapse when the 
presynaptic voltage exceeds Vth ■ We assume the synapses 
are fast, i.e., the time constants of the synapses are taken 
much less than period of the spiking neuronsSS. In this 
regime the coupling terms can be approximated by pul- 
satile currents (see below). 

Given the choice of model parameters in TABLE. HI an 
uncoupled WB neuron exhibits periodic spiking with a 
period of about 16ms for li = IfiA/cm?, and a HH neu- 
ron fires with a period of about 14ms for li = lOfiA/crr?. 
Throughout this paper we have taken this typical values 
for I2 and then inhomogeneity is imposed by Ii — I2 + AI 
with A/ > 0. consequently, wi = a;2 + Aw and the first 
neuron has a larger natural firing rate with Auj > 0. 

Dynamics of limit cycle oscillators can be approxi- 
mated by the equation describing the evolution of the 
averaged phases. For an analytic inspection we use this 
approximation to quantify the evolution of the membrane 
voltage of two regularly firing neurons connected by pulse 
like couplings with Winfree type oscillators'^^: 

e^=UJ^+ 9^JQm E„ ^{t " " ^j) , (2.4) 
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where 5 is the Dirac's delta function describing the pul- 
satile interaction, and gij is the strength of the synapse 
which can take both positive and negative values to de- 
scribe excitatory and inhibitory synapses, respectively. 
The phase reset curve Q{0) determines the response of 
the neuron to the incoming pulse. Delay time between 
the firing of the neuron j and elicitation of a postsynap- 
tic pulse in neuron i is denoted by Tij. uJi is the natural 
firing rate of the neuron (when isolated). 

We have mainly studied the synchronization of the neu- 
rons by recording intervals between spiking of two neu- 
rons. Furthermore we have measured synchrony by cal- 
culating the zero lag cross-correlation between phases of 
two neurons: 



Pl2 = 



jo imit)) - {Oi{t)){0iit)) 

^Var{9i)^Var{e2) ' 



(2.5) 



Where Var is the variance. For the limit cycle oscillators 
of WB and HH type, interval between two successive ac- 
tion potentials defines a complete cycle, and the phase 
increase during this time amounts to 27r. Then we can 
assign a value to the phase by linear interpolation be- 
tween two action potentials: 



27r 



(2.6) 



where U is the time of the last spike and T is the inter- 
spike intervalii. Advantage of using the correlation func- 
tion as a measure of phase locking is that it can discrim- 
inate between uncorrelated firing and antiphase locking 
since it assigns a zero value to the former and negative 
values to the later case. Order parameters of the Ku- 
ramoto type^^ used by other authors (see e.g. [20]) take 
zero value for both the antiphase and uncorrelated fir- 
ings. 

The differential equations were solved using a fourth- 
order Runge-Kutta method3§. with the step size dt = 0.05 
of real time. Each simulation typically lasts 200000 iter- 
ations (lOsec). 

Phase reset curves: Phase reset curves determine 
the responses of neurons to brief stimulations. As illus- 
trated in Fig. [TJ a pulse of duration 1ms and strength 
h in different phases 9 is imposed after the ith spike, 
and the timing of the next spike U+i is recorded. The 
normalized phase reset Q{0) is defined as 



Q{0) = r(i 



t+i 



u 



T 



(2.7) 



where the phase is defined by Eq. 12.61 and h is the 
strength of the pulse which can be positive or negative 
for excitatory and inhibitory pulses, respectively. 

In Fig. [21 the phase-response curve for positive (ex- 
citatory) and negative (inhibitory) pulses are shown for 
the WB and HH neurons. The most notable difference 
between the two models is that PRC is strictly posi- 
tive for the WB model, whereas it exhibits both positive 
and negative regions in the HH case^^. Hence, a small 



9 = 



6 = 271 




time 



FIG. 1: Illustration of the measurement of the phase response 
curve of the model neuron. The solid curve indicates the 
voltage of the unperturbed neuron, while the dashed curve 
indicates the voltage after a brief perturbation is applied at 
time t. 
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FIG. 2; Normalized phase Reset Curves for the WB model 
(top) as a type-I neuron and the HH model (bottom) as 
a type-II neuron are plotted for excitatory and inhibitory 
pulses. Horizontal axis shows the phase at which a pulse 
of duration 1ms and the amplitude 0.05 is imposed on the 
model neuron. 



depolarizing perturbation always results in an advance 
in firing for the WB model, whereas it may either ad- 
vance or delay spiking in the HH model, depending on 
when exactly the perturbation is delivered. Both the WB 
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FIG. 3: A schematic of the map which is used for analytic 
investigation. The evolution of the averaged phase of two 
neurons is shown by sohd and dashed hnes, respectively. T 
is period of firing in phase locked state and nj is the delay 
time. Dynamical variable of the map is the phase difference 
between nth spikes of two neurons depicted by A0„. 

and HH models show maximal sensitivity to the exter- 
nal perturbation far away from the spiking events. In 
addition, the HH model shows more evidently that right 
after the spike, the system is rather unresponsive to in- 
coming perturbations, due to refractoriness. The results 
for inhibitory and excitatory pulses do not show any qual- 
itative difference. 



locked state: 

TAu} = g2iQioJ2T2i - AO) - guQiu^iT^ + AO). (3.3) 

Note that in general the network period T is not con- 
stant and depends on the other parameters. Intuitively, 
this equation expresses that the phase difference arose 
from the mismatch over a period should be balanced by 
the mutual excitations (or inhibition) for phase locked 
activity. For given values of parameters, a solution of 
this equation, if exists, determines the phase lag of fir- 
ing of two neurons in locked state. An upper limit for 
the mismatch which the system can tolerate in a phase 
locked solution is given by 

TAuj < Max[g2iQ] - Min[gi2Q]. (3.4) 

Stability condition: Assuming a small perturbation 
on the phase lag in locked state A6'„ = Ad + Cn , we get 
the linearized map 

Cn+l = C«[l + gi2Q'{uJlTi2 + AO) + g2lQ'{i^2T21 " AO)], 

(3.5) 

where Q'iO) — — — . Stability condition is 
dO 

gi2Q'{uJiTi2 + Ad) + g2lQ'{uj2T21 - Ad) < 0, (3.6) 

which guaranties |Cn+i| < |Cn|- 



III. RESULTS 



A. Analytic results 



Out of the causal limit (when the incoming pulses do 
not elicit action potentials instantaneously), we construct 
an iterative map by sampling the phase of the Winfree 
oscillators, defined in Eq. 12. 4[ at the time of firing of the 
faster spiking neuron, neuron 1. Assuming Oiit") = 6i{n) 
we have 



L(n + 1) = 6*1 (n) + ujiT + gi2Q{uJiTi2 + A6i„), 

2{n + 1) = 62(71) + UJ2T + g2lQ{0J2T2l - AOn). 



(3.1) 



Here A6'„ is the phase difference of the neurons at the 
instance of nth spike of the high frequency neuron (see 
Fig. in]), and T is the period of the firing of the high 
frequency neuron which in general is different from its 
natural period of the firing. In the locked state T is the 
period of the firing of the both neurons in the network 
(see below). 

A I : 1 locked state is characterized by a fixed point of 
the map 



AO. 



n+l 



AOr, + TAuj 



ffi2Q(wiri2 + Adn) 

g2lQ{L02T2l- AOn), (3.2) 



which gives an implicit equation for the phase lag AO in 



1. Special Cases: 

a- Symmetric configuration: For the bidirectional 
symmetric couplings 512 = 321 = g and T12 — T21 = r, 
and identical neurons uji ~ uj2 ~ (jJ Eq. 13.31 gives 



Q{u:t + AO) = Q{ujT- A9). 



(3.7) 



This equation has always a synchronous solution ~ 
regardless of the functional form of Q. But since Q is 
a periodic function of phase, assuming it is continuous, 
other solution also exist which depends on the functional 
form of PRC. Stability condition for synchronous (in- 
phase) state is gQ'{u)T) < consistent with previous 
results^' 

b- Symmetric couplings with dissimilar neu- 
rons: For the bidirectional symmetric couplings gi2 ~ 
321 = g and ri2 = r2i = r, Eq. 13.31 gives 

TAuj = g[Q{uj2T - AO) - Q(wit + AO)]. 

Assuming small mismatch, we can linearize Q and deduce 
an explicit equation for small phase lag: 



AO 



-Auj\ 



T 



'2 2gQ'{u}T)' 



(3.8) 



For zero mismatch, synchronous solution exists indepen- 
dent of the delay and synaptic strengths. For nonzero 
mismatch, inphase firing is only possible if the terms in 
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FIG. 4; The phase difference between the spikes of two neu- 
ron is plotted against the delay in communication for type I 
and type II model neurons. The model neurons are iden- 
tical and are of the Winfree type introduced by Eq. 13.11 
with iJi — UJ2 — 1. The synapses are excitatory with 
912 = 321 ~ 0.5. The slope of phase reset curves are also 
shown by shaded regions, green regions show the ranges where 
synchronous firing is stable gQ'{u}T) < 0. 



brackets cancel each other (note that second term is nega- 
tive due to stability condition given below). Yet for small 
mismatch, near isochronous firing is possible for a range 
of delay times since the phase lag changes proportional 
to delay with the rate Auj <^ 1. Stability condition given 
in Eq. 13.61 for the symmetric connections gives: 



glQ'icJiT + AO) + Q'iuj2T - AO)] < 0. 



(3.9) 



F4or small mismatch, if a synchronous state exists, it is 
stable if gQ'iujT) < 0. 

c- Near symmetric configuration: In a near sym- 
metric configuration when the neurons are not identical 
and the connections are not exactly symmetric, such that 
Auj = oji —102, Ag — g2i — gi2 and At — T21 — T12 are 
small parameters, the phase lag in a near inphase state 
can be given by a linearized approximation: 



A6 



1 



2gQ'' 



-{T + gTQ')Auj - QAg - gojQ'Ar], (3.10) 



where Q — Q{ujt) and Q' = Q'{u!t). Together with the 
stability condition gQ' {ujt) < this equation determines 
deviation from synchrony for given amounts of hetero- 
geneity in the neuronal and synaptic parameters. This 
equation shows that in principle it is possible to tune the 
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FIG. 5: Same as the Fig. |4] the phase difference between 
the spikes of two neuron is plotted against delay time, here 
for two nonidentical neurons with a mismatch shown in the 
plots. Lower plot shows the phase difference in respect to the 
mismatch for a sample value of delay time. The red solid line 
in the lower plot is the result of the analytic approximation 
Eq. 



phase lag by the changes in variable parameters, namely 
synaptic strength. 

d- Unidirectional coupling: Assuming a unidirec- 
tional coupling gi2 — necessary condition for phase 
locking is TAuj < Max{g2iQ)- With excitatory synapse 
.921 > 0, this equation states that over each period, phase 
shift exerted by excitatory synapse on the low frequency 
neuron should be strong enough to compensate for the 
difference in natural frequency. For inhibitory synapse 
g2iQ should have a positive lobe for existence of phase 
locked state, i.e., just type II neurons can be locked to- 
gether. The phase lag in the locked state for unidirec- 
tionally coupled neurons can be explicitly found 



Ae = LJ2T21 - Q-\ 



TAlo^ 

521 



(3.11) 



As can be seen phase lag changes linearly with r and if 



6 



Aw 7^ 0, can be changed by i.e. in this case plastic- 
ity can tune the phase difference. On the other hand in 
absence of mismatch Aw = 0, phase lag is independent of 
synaptic strength and is solely determined by firing rate 
and the delay time. 



B. Numeric results 

We first give numeric results for two pulse coupled 
Winfree oscillators to check the analytic results given 
above. Phase reset curves for type I and type II oscil- 
lators are assumed as 1 — cos{9) and sin(9), respectively. 
In Figs, m and [5] we have shown the time lag between 
the spikes of two bidirectionally coupled model neurons 
in symmetric and near symmetric configurations. For a 
fully symmetric case, when both the neurons are identi- 
cal and couplings are symmetric, the results match the 
results of other authorai^. With the excitatory synapses 
the neurons fire in synchrony (with zero phase lag) for 
the range of the delay time in which the PRC has a neg- 
ative slope. This is consistent with Eq. I3.3l where A9 = 
is a trivial solution regardless of other parameters, along 
with the stability criterion gQ'{LdT) < 0. 

Any deviation from the symmetric case makes the 
phase locked state to change with other parameters as 
delay, synaptic constant and firing rate. For example 
with nonidentical neurons modeled with different firing 
frequencies, Eq. 13.31 has not trivial solutions and the 
phase lag of firings in locked state changes with delay 
(Fig. [5]). For small mismatch a linear approximation can 
show the dependence of the phase lag on other parame- 
ters as is given in Eq. 13.81 Numeric results shown in Fig. 
[5] conform with analytic results for small mismatches. It 
can also be seen that transition from the inphase state 
to antiphase, is mediated by an unlocked region in which 
the neurons fire independently. This unlocked region in 
centered around the point in which derivative of PRC 
changes sign and neither inphase and antiphse solutions 
are stable. For the fully symmetric configuration it oc- 
curs just in the point of zero slope of PRC, but the re- 
gion grows with increasing mismatch. In another point of 
view, the results show that the range of mismatch which 
the system can tolerate and remain in a phase locked 
state, depends on the slope of phase reset curve. Hence 
the tolerance of the locked state to the mismatch depends 
also on delay time: Near the regions in which slope of 
PRC changes sign, the system is vulnerable to the het- 
erogeneity and nonidentical neurons can fire neither in 
synchrony nor in phase locked manner. 

We have repeated numeric experiments for conduc- 
tance based neuronal models, WB and HH as type-I and 
type-II neurons, respectively. Results shown in Figs. [S] 
and [7] qualitatively match the results for Winfree oscil- 
lators. Synchronous state is not structurally stable and 
phase lag between the firings changes proportional to the 
mismatch for small mismatches (see the middle plot in 
Fig. [5]). Wang et al. have shown that with large delays 
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FIG. 6: Phase difference between the spikes of two neuron 
is plotted against delay time, for two bidirectionally coupled 
WB neurons with excitatory (top) and inhibitory synapses 
(bottom). The results are shown for both identical neurons 
and in presence of mismatch. Input currents are /2 = 1 and 
Ji = 1 + A/. The currents are chosen such that natural 
period of firing for the slower neuron is about 15ms. Middle 
plot shows the phase difference vs. mismatch for WB neuron 
for a sample value of delay time. 



synchronous firing is more stable against mismatch, com- 
paring to small delay timesi^. Our results for HH model 
(Fig. [7]) are consistent with this result but in a general 
point of view dependence of phase lag to the mismatch 
is determined not only by delay time but also by other 
parameters as coupling constant and the slope of PRC 
for a given delay time (see Eq. 13. 8p . So the form of PRC 
also determines the sensitivity of the synchronous state 
to the mismatch, for example, for a canonical form of 
PRC of type-II which is symmetric around 9 = n, phase 
lag changes similarly for small and large delays (see Fig. 

ED. 

In the formalism presented in previous section, both 
PRC and its derivative always appear as their multipli- 
cations with corresponding coupling constant, i.e. as gQ 
and gQ' ■ If the PRCs maintain they form for inhibitory 
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FIG. 7: Phase difference between the spikes of two neuron is 
plotted against delay time, for two bidirectionally coupled HH 
neuron with excitatory (top) and inhibitory synapses (bot- 
tom). The results are again shown for both identical neurons 
and in presence of mismatch. Input currents are 72 = 10 for 
WB neuron and = 10 for HH neuron and Ji = 1 + AJ. 
The currents are chosen such that natural period of firing for 
the slower neuron is about 15ms. Note that with the applied 
mismatch in input currents shown in this figure and Fig. [G] 
relative mismatch is equal for both models A/// = 0.05. 



pulses 5 < 0, with the transformation g —g the gen- 
eral result for existence of phase locked solutions (Eq. 
13. 3p remains invariant with the simultaneous permuta- 
tion transform 10 2. But the right hand side of the 
stability condition (Eq. 13.61) for inphase and antiphase 
solutions changes sign and therefore the regions of sta- 
ble synchronous solutions loose stability and vice versa. 
For identical oscillators, for example, this means that the 
regions of inphase and antiphase solutions interchange. 
Results shown in Figs. [5] and [7] support this results: In 
both the figures the regions for inphase firing (near in- 
phase firing for nonidentical neurons) are almost inter- 
changed with regions of antiphase firing (near antiphase 
firing for nonidentical neurons) . Specifically, with excita- 
tory synapses and small values of delay time, typical type 
I neurons fire in antiphase manner while type II neurons 
fire synchronously. 

According to Eq. 13.101 it is possible to maintain syn- 
chrony, if different source of inhomogeneity show com- 
pensating effects. We have tested this hypothesis by 
considering nonidentical oscillators coupled by unequal 
strength synapses. Figure [5] shows that for a given delay 
time if synchrony is stable for symmetric configuration. 
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FIG. 8: Multiple sources of inhomogeneity can compensate for 
each other and make synchrony. Top, A Density plot of cross 
correlation coefficient is shown against mismatch in firing rate 
and mismatch in coupling constants for two type-II Winfree 
oscillators. A temperature-map color scheme is used to assign 
warmer colors to higher values of correlation which indicate 
inphase synchrony. Delay is chosen such that the symmetric 
configuration shows inphase synchrony. It can be seen that 
nonidentical neurons can be synchronized by unequal coupling 
constants. The thick line shows result of Eq. IS.lOl for inphase 
synchrony. Bottom, the result presented for two HH neurons 
shows similar qualitative behavior. 



simultaneous changing of asymmetry parameters (synap- 
tic strengths and firing rates) can maintain synchrony. 
Note that for the given example, the synapse from the 
high frequency neuron to low frequency neuron should be 
stronger to bring the neurons back to synchrony; an in- 
tuitively reasonable result. For the HH neurons, interest- 
ingly, the region of near inphase firing is larger comparing 
to Winfree oscillators but yet the results hold qualita- 
tively. 
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IV. DISCUSSION 

In most studies on the synchronization of two de- 
layed coupled neurons, a symmetric configuration is 
assume d^'^'^^ . This means that the neurons are assumed 
identical and the couplings are considered with exactly 
equal weights and equal delay times. To step beyond this 
seemingly unrealistic assumption, in this study we have 
considered an almost symmetric configuration where ei- 
ther the neurons (nodes) or the coupling terms (edges) 
have slightly different parameters. Our analytic results 
for two pulse coupled phase oscillators show that any 
small inhomogeneity in a parameter, when imposed on 
the system, can change the synchronous firing to an al- 
most synchronous firing where the neurons fire with a 
small phase lag which is proportional to the parame- 
ter of inhomogeneity. These parameters in our minimal 
model are the natural firing rate of the neurons, and the 
strength and the delay time of two directed couplings. If 
for a symmetric system synchronous firing is stable, then 
any small asymmetric deviation in firing rates, synap- 
tic strengths and delay times exerts a phase lag between 
firing of two neurons. 

As is evident from Eq. 13.101 the functional form of 
phase reset curve of the neurons and the strength and 
the type of synaptic coupling (inhibitory /excitatory) de- 
termine deviation from inphase firing. As the main out- 
come of this study we have shown that different sources 
of inhomogeneity when are present in a system together, 
can cancel each other and bring the system closer to syn- 
chrony. For example, for nonidentical neurons, it is re- 
cently shown (see Ref. [19]) that synchronization is only 
possible in the causal limit (when the excitations elicit 
an immediate action potential in postsynaptic neuron) 
with equal delays. In this case synchronization is not 
possible with short delay times and a symmetric increase 
in synaptic conductances is not a effective strategy to 
decrease the time lag. Instead, we have shown that an 
asymmetric change in synaptic strengths can bring the 
neurons back in synchrony (when the two first terms in 
the brackets in the right hand side of Eq. 13.101 cancel 
each other). Interestingly, spike timing-dependent plas- 
ticity (STDP) exerts an asymmetric change in bidirec- 
tional synaptic couplings, but with the classical profile, 
STDP breaks neuronal loop and leads to the divergence 
of mutual couplings such that they are just restrained 
by limiting values^. The outcome of operation of STDP 
in this case will always be a unidirectional coupling and 
hence, STDP with classical profile is not a good candidate 
to compensate the asymmetry in the system. It needs 
further investigation to check if nonlinear profiles of spike 
timing-dependent plasticity can show such self-regulatory 
effect on the inhomogeneous systems'^'^ . Interestingly, ax- 
onal conduction velocity can also be modulated^. It is 
not clear if there is a self-tuning mechanism to change 
the conduction velocities but in principle variable delays 
also can role as another regulatory parameter to balance 
the system and bring the neurons closer to synchronyi^. 



Introducing a relay neuron in the midway between two 
distant neuron is a sensible proposal to extend the do- 
main of stable synchrony in the parameter space22i2i. 
Indeed in this case both antiphase and inphase firing of 
adjacent neurons, lead to simultaneous firing of wing neu- 
rons and so synchrony will be stable for almost all values 
of delay time. But again, in the full parameter space, 
this synchrony can be achieved in symmetry manifold on 
which all the parameters of two wings are equal. It is 
the subject of future study to investigate how deviation 
from symmetric configuration affects the synchronization 
properties of such a system. 
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Appendix A: The model neurons 

Wang-Buzsaki (WB) Model: The steady-state ac- 
tivation moo and the rate equation for the inactivation 
variable h in the expression for sodium current and rate 
equation for the activation variable n in the expression 
for potassium current are given respectively as follow: 



m = rriooiy) = 
dh 



— = (j,[ah{V){l-h)-Pn{V)h], 



(Al) 



n 
dt 



{a„iV){l - n) - MV)n]. 
The rate constants for TOqcj h, and n are: 



am{V) = -0.1{V + 35)/(exp(-0.1(y -I- 35)) - 1), 
/3„(y) =4exp(-(T/-H60)/18), 

ahiV) = 0.07cxp(-(y + 58))/20), 
Ph{V) = l/(exp(-0.1(y + 28) + 1), 

an{V) ^ -0.01(1/ + 34)/(exp(-0.1(y + 34)) - 1), 
MV) ^ 0.125exp(-(F-|-44)/80). 

(A2) 

Ho dgkin- Huxley (HH) model: The rate equations 
for the activation variable m and inactivation variable h 
of the sodium expression, and n, activation variable of 
potassium, obey the differential equations: 



dfn 

— = [a„,{V){l ~ m) ~ l3,„{V)ml 
dt 



dt 

dr. 

lit 



= [ah{V){l - h) - Ph{V)h], 
[a„(T/)(l-n)-/3„(l/)n]. 



(A3) 
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The rate constants for m, h, and n are: 





= (2.5-0.iy)/(exp(2.5 


-O.IV) - 1) 




= 4 exp [—V/io), 






= 0.07 exp (-y/20), 




^h{V) ■■ 


= l/(exp(3-0.iy) + l), 






= (0.1-0.0iy)/(exp(l- 


o.iy)-i), 




= 0.125 exp (-y/80). 
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